The technology used for data acquisition will affect how 3D data should be processed prior to the morphometric analysis. Digitizers directly collect coordinates of points from the physical object. Hence, the data retrieved is usually ready for analyses by GM toolsets after little or no processing. As a counterpart, errors or modifications of the landmark configuration often imply to recollect data and hence to access the specimen again.

On the contrary, data registered from surface and tomographic scanners require tedious software postprocessing to obtain a 3D digital model of the specimen and collect morphometric data on it. High-resolution scans yield large amounts of data resulting in very dense point clouds. They require important decimation before export to morphometric software. Carefully selecting appropriate cleaning, alignment, merging and decimation options is critical in order to preserve the geometrical features of the original object in its digital model, while at the same time significantly reducing the memory load. Various software packages offer solutions to collect geometric morphometric data on the digital mesh model. Software for mesh processing or data visualisation (VSG Avizo Fire 6.1, http://www.vsg3d.com/; Geomagic Studio 11, http://www.geomagic.com/) provide tools to register point coordinates by hand or semi-automatically. A limited number of software packages, however, are explicitly oriented for 3D geometric morphometric data acquisition and/or processing and provide varying levels of user-friendliness (Edgewarp3D 3.31–Bookstein and Green, 2002; IDAV Landmark Editor–Wiley et al., 2005). Though scanning is time demanding, digital models have the advantage to allow for easy and fast correction or management of datasets as the morphometric acquisition step can be saved in a file, to be reopened and edited.

At this point, we highlight that all these methods result in a representation of the structure(s) of interest made of discrete 3D landmark coordinates, likely to be used in quantitative analyses of shape variation.

The ‘template’ is a landmark configuration used as the reference for any new specimen added to the study. It is composed of every landmark used in the analysis with its associated type–true landmark and either curve or surface semi-landmark. The choice of the type, number and relative location of landmarks is hence a critical step in creating a template configuration, that is highly influenced by the aim of the study and the nature of the sample.

Homologous (“true”) landmarks define a direct correspondence between parts of different objects: homology is in this case not the biological homology, but a geometrical one (sensu Jardine, 1969). True homologous landmarks belong to the types I and II of Bookstein (1991) that correspond to juxtapositions of tissues and maxima of curvatures respectively. By opposition, type III corresponds to weakly defined points (the more distant from a reference by instance) and to outline points. Both are regrouped into the name of semi-landmark that recalls that they are defined only in one direction and that some of their coordinates are deficient (Bookstein, 1997 and Gunz et al., 2005).

Semi-landmarks demand different approaches than true landmarks to be used in shape comparisons. Two strategies can be employed.

In the first procedure, introduced by Bookstein (1997) (see also Gunz et al., 2005), non-homologous landmarks can be transformed into operational landmarks using geometrical and mathematical criteria. This is achieved by allowing them to slide, whether along predefined curves or on surfaces, while optimizing a criterion such as minimizing the bending energy parameter used by the thin plate splines (Bookstein, 1997, Gunz et al., 2005 and Mitteroecker and Gunz, 2009). True homologous landmarks are mixed with non-homologous ones in order to introduce spatial constraints (Gunz et al., 2005). At the end of the process, all landmarks, whatever their types, are superimposed and treated in the same way true homologous landmarks are. To perform this operation, curve and surface semi-landmarks must be transferred from the template configuration onto the target specimen by ‘warping’. In this procedure, the new target specimen is defined by its true landmarks coordinates. Based on this different configuration of true landmarks, the thin plate spline is a one-to-one spatial interpolation function that warps semi-landmarks into their new position. Once warped, semi-landmarks must be projected onto their respective target curves and surfaces to proceed to spline relaxation. Finally, as explained in detail in Gunz et al. (2005), semi-landmark sliding must be repeated iteratively using each newly formed consensus until convergence is attained.

A second approach uses mathematical functions to directly describe the whole set of points. This has been done with 2D and 3D outline points (see Fig. 1 and Fig. 2 for examples of 3D outlines) by using, for example, Fourier decompositions (Baylac and Friess, 2005, Ferson et al., 1985 and Kuhl and Giardina, 1982) or eigenshape analysis (McLeod, 1999). At one point, steps are included that allow for some form of registration, including size normalisation, either of the raw coordinates or directly applied to the function parameters (Kuhl and Giardina, 1982). At the end of the process, the parameters of the functions are used as shape variables in the statistical analyses. Invertible functions have been favoured since their parameters may be used to reconstruct the shapes allowing for straightforward visualisations of the shapes.

For homologous and sliding landmarks, a single morphometric operation, the Procrustes superimposition (Rohlf and Slice, 1990) is used to provide a common referential system that makes point coordinates comparable. This step is necessary because, unlike distances, point coordinates are not invariant to the simple geometric operations of translations and rotations: measured distances remain identical if one translates and rotates the objects, whereas point coordinates do not. The superimposition procedure will eliminate from the original coordinates the components related to the translations and rotations, adding a further size normalisation operation. Mathematically, this superimposition operation minimises the overall sum of squares between all homologous landmarks of a set of objects. Geometrically, it is a translation and a solid rotation of each object toward a reference object, usually iteratively calculated as the mean object or consensus. The result will be a set of superimposed landmark coordinates that have the same size and that summarize the geometrical differences between objects. The superimposition process therefore realizes the extraction of size and shape (here the superimposed coordinates) terms from the overall form of an object, following Needham's (1950) equation: form = size + shape. These size and shape terms may be simply visualized but more interestingly can be subjected to multivariate statistical analyses (Dryden and Mardia, 1998 and Krzanowski, 1988).

Size normalisation is done by dividing each set of point coordinates of an object by its size. Size is estimated by the centroid size (Bookstein, 1991), the square root of the sum of squared distances between the centroid of an object and its landmarks. Superimposed coordinates are therefore ratios and size normalisation corresponds to a simple scaling step, that leads to objects of uniform size but does not modify their individual proportions. Nevertheless, biologically, size and shapes are rarely independent. Moreover, their relationships define one important component of the biological variability, namely the allometries (see Klingenberg, 1996). Readers interested by the important conceptual and practical consequences of the use of ratios are referred to Darroch and Mosimann (1985), Mosimann and James (1979) and Bookstein (1996) together with Small (1996) and Dryden and Mardia (1998), which introduce the complexity of shape spaces (see also Rohlf, 2000 for useful and visually illuminating examples).

For some reasons, detailed in Gunz et al. (2005, pp. 86–88), statistical inferences using semi-landmarks, whether of the surface or of the ridge type, are not a simple straightforward extrapolation of the usual inferences done using the Procrustes coordinates. Nevertheless, exploratory analyses like Principal Components Analyses (PCA) or Partial-Least-Squares (PLS) are always feasible. In this article, the description of the shape variability uses exclusively PCA.

The value of GM does not only reside in its great statistical power. It also resides in the simple yet powerful visual interpretations it provides. Visualisations may be simply direct comparisons of superimposed shapes. More generally, the analyses of shape variability are frequently done in relation to predefined factors (size in the case of allometry, groups of objects, archaeological beds, palaeontological strata, ecological or environmental parameters) that may be related to particular directions of variability into the shape space. Alternatively statistical directions like factorial axes, principal components or canonical axes may be also specifically calculated from the shape covariance matrices. Visualisations may also be done along the directions of maximum variances between one set of shape parameters and one set of variables of interest using 2-Blocks Partial Least-Squares (2B-PLS; Rohlf and Corti, 2000). This permits analysis and visualization of shape changes associated with a set of variables characterising past or present environments, ecotones, ecological niches, functional or physiological parameters (Adams and Rohlf, 2000).

All these visualizations of either 2D or 3D landmarks and semi-landmarks estimate the statistical relationships between shape and the specific directions of interest (DOI) using multivariate regression (Dryden and Mardia 1998; Monteiro, 1999). By this methodology, the particular landmark configuration that corresponds to a single value picked along the DOI may be calculated and visualized. Usually both extremities along the DOI are contrasted to illustrate the overall shape differences (see “Results”). Multivariate regression is useful whatever the nature of the variables is: landmarks, semi-landmarks, outline or surface function parameters. In the case of a shape described by parameters of mathematical functions, the inverse function allows for the calculation and visualization of the coordinates of reconstructed shapes.

The first dataset deals with large scale evolution and locomotory adaptations among extinct and extant theropod dinosaurs (bipedal, generally carnivorous; including birds). The history of this diversified group dates back to Carnian (Late Triassic > 230 My) and covers major geological and evolutionary events. Bipedalism was present since the origin of Theropoda (and perhaps other Dinosauria), while the > 9500 species of modern birds today remain the sole obligatory, striding bipeds with the exception of humans. Moreover, their terrestrial locomotor apparatus maintained an overall similar organisation, though it underwent a gradual sequence of modification that has motivated numerous studies of anatomy and function (Gatesy, 1995, Hertel and Campbell, 2007, Hutchinson, 2001, Hutchinson, 2004a, Hutchinson, 2004b, Hutchinson and Gatesy, 2000, Rasskin-Gutman and Buscaglioni, 2001, Reilly, 2000 and Rubenson et al., 2007).

The aim of the broader study is to analyse shape variation of the pelvic girdle with a focus on the geometry of hip joint and antitrochanter (articular surface adjacent to the acetabulum, Fig. 1). To investigate how hip form and function evolved, the sample must be representative of the morphological, ecological and taxonomical diversity found among Theropoda. Hence, this dataset is remarkable because of the important variability it comprises. As a result, finding homologous landmark locations can be problematic and no single specimen can accurately represent the morphology of all others to create a template configuration (Fig. 1).

To illustrate such analysis, the Theropod Dataset presented here comprises 3D models of five pelvic girdles (os coxae). Extant animals are represented by four randomly selected genera of modern birds belonging to three different clades: Coturnix (Phasianidae, Galliformes), Calonetta (Anatidae, Anseriformes), Troglodytes (Trogloditidae, Passeriformes) and Taeniopygia (Estrildidae, Passeriformes). These genera range from very small (9 cm long, 10 g) to small (20 cm long, 180 g) birds. Models were obtained using X-ray microtomography, with isometric voxel resolutions ranging from 22 μm to 64 μm. Segmentation was performed using Avizo (VSG Avizo Fire 6.1).

Allosaurus fragilis is the only representative of non-avian theropods included in this study. Scanning was performed at the Museum of the Rockies on fossil specimen MOR 693 (‘Big Al’) with a Modelmaker scanner and postprocessed with Paraform Viewer.

We placed a total of 936 landmarks on each specimen (Fig. 1): 44 anatomical landmarks, 280 sliding semi-landmarks on 3D curves and 612 sliding semi-landmarks on 3D surfaces. Anatomical landmarks were placed on homologous locations delimiting the overall shape of the pelvis, as well as circumscribing the endpoints of curves. The curves represent clearly defined ridges, either delimiting the extent of the preacetabular ilium or the articular surface of the acetabulum and antitrochanter. Surface sliding semi-landmarks form a dense mesh over the articular facets of both hip joints.

The second dataset deals with skull evolution in populations of shrews. One of the most noticeable particularities is the homogeneity of their global form across time and geography (Churchfield, 1990). Among shrews, Crocidura russula is morphologically wellconserved throughout its geographical range. However, previous studies using 2D GM on French populations showed tiny but significant morphological differences related to geography in general and insularity in particular.

The aim of the current study is to characterise and quantify more precisely morphological variations on the whole skull (Fig. 2) using 3D surface GM in order to explain these variations linked to genetics or/and environmental factors.

This dataset is remarkable because different specimens are very similar and potential differences very tiny. Nevertheless 3D surface GM allow us to describe both globally the form of the skull and specifically the different regions of the skull involved in sensory functions like the ocular, the auditory, the nostril or the masticatory regions. Another consequence is that one specimen can potentially be used to create a template configuration efficient for the whole dataset.

To give an example of such application, the Shrew Dataset is composed of five specimens from five different populations. These specimens have been scanned using X-ray microtomography with a resolution of 20 μm. As for the precedent dataset, data have been segmented with Avizo (VSG Avizo Fire 6.1).

996 landmarks have been used to describe the skull. Most of them are spread over the dorsal side (Fig. 2). They consist of 39 anatomical landmarks, 157 sliding semi-landmarks on 3D curves and 800 sliding semi-landmarks on 3D surfaces. Most anatomical landmarks correspond to intersections of different bones or teeth and can be considered to belong to type “1” homologous landmarks, according to Bookstein (1991). Curves are located at the periphery of large bones like the parietal, maxillar and nasal. Sliding semi-landmarks on 3D surfaces were placed over the dorsal part of the skull from the nasal to the parietal bones, where no true landmarks could be found. Semi-landmarks were also used to describe the mandibular fossa region which, as most articular surfaces is devoid of any true landmarks.

All these choices have been done regarding two aims, first to describe the global form of the skull by a relatively dense covering of semi-landmarks, and second to deal with regional part of interest of the skull, most of them corresponding to sensitive regions.

As the Generalised Procrustes Superimposition (GPA) relies on a landmark-based discrete representation of biological shape, the design of an appropriate template is a crucial step in the workflow. It is pointed out by Gunz et al. (2005) that this definition of the relative location and repartition will dominate the resulting statistical processing and visualisation of the data.

In their paper, Gunz et al. (2005) explain how they proceed to obtain the landmark configuration that is used as the template. To do so, they use a randomly chosen specimen of their sample as reference. True landmarks and curve semi-landmarks are digitized directly onto it, while surface semi-landmarks are obtained from a dense point cloud digitized on the specimen. Decimation–controlled reduction of the number–of the redundant points of this cloud produces a more sparse and relatively regular pattern that can be used as surface semi-landmarks. As they specify, similar results can be obtained by using the dense mesh resulting from surface or volume scanning.

In our case, the different nature of the two dataset requires different approaches for the design of each template configuration. This leads to the introduction of a novel method to design templates for sliding semi-landmark GM analysis, adaptable to the level of shape variability in the sample. During the warping of the template onto the target specimen, as the configuration of true landmarks interacts with semi-landmarks, much of the original geometry of the template is retained for curve and surface semi-landmarks. As a result, for semi-landmarks located on morphological features geometrically similar on the template and the target, warping will efficiently bring them close to the homologous target area. This proves particularly useful in the case of samples with very homogeneous specimen morphologies, as little to no hand processing has to be performed after warping, before to project and slide the semi-landmarks. On the contrary, for features with significantly different shape on both objects, warping may leave the semi-landmarks relatively far from the target area they should be projected on. These cases require difficult and time-consuming hand processing to correct the position of offset semi-landmarks and attain a satisfying match onto the target specimen, so that the semi-landmarks can be correctly projected and relaxed.

As previously noted, contrarily to the Shrew Dataset, no single specimen from the Theropod Dataset displays morphological properties that can efficiently be warped into the geometry of every other specimen. Hence, we designed both templates using polygonal modelling tools in Maya (Autodesk^® Maya 2010) to create each template with a different level of geometrical complexity. With this method, instead of designing the Theropod Dataset template on one random specimen, we were able to create it with a much lower level of complexity (Fig. 3). It was modelled using simple 3D geometrical primitives–planes, cylinder–and controlled modelling tools, to produce a basic idealised shape with equidistant curve semi-landmarks and grid-like surface patches. Because the template shape is less complex, warping the semi-landmarks retains simpler geometrical features. As a direct benefit, fewer and faster manual adjustments are required before the warped semi-landmarks closely match their new target position.

On the other hand, the Shrew Dataset presents a very homogeneous form. Hence, we could model a reference mesh with greater geometrical complexity based on the morphology of a specimen of the dataset (Fig. 3). As a result, we obtained a template that can easily be warped and projected onto other target specimens, in a similar way to what is proposed in Gunz et al. (2005). However, polygonal modelling has other benefits as it provides direct control over the number and the individual location of semi-landmark points. Moreover, the simple tools used to add or refine divisions and vertices allow one to create very regular tessellations with evenly spaced points covers. This modelling of the reference mesh into a regular pattern permits an efficient management of local point density (Fig. 3) and to subdivide the template into regions of interest for further analysis. The differences in the warping step resulting from the distinct approach used for each dataset are illustrated in Fig. 4.

We use Edgewarp3D 3.31 (Bookstein and Green, 2002) to warp our reference configuration onto each specimen and perform iterative semi-landmark relaxation to their final target positions. Fig. 4 illustrates this step for both datasets. After an initial Procrustes superimposition and computation of a consensus configuration, semi-landmark sliding is iteratively repeated using the new computed consensus as template until convergence is obtained (Gunz et al., 2005). The resulting data is exported to perform GM in the graphical and statistical R environment (http://cran.r-project.org/) using the Rmorph library (MB). Each dataset is then superimposed using a Generalized Procrustes Analysis. Objects are made symmetrical using the procedure for object symmetry of Mardia et al. (2000). Superimposed coordinates are subjected to an unscaled Principal Component Analysis (PCA). Visualization are done using multivariate regression and links between selected points are defined in order to help visualize the objects. The results are used to produce warped surfaces and visualisations using two programs, namely SurToPly and ColorMorphoPly, which are still in ongoing development by one of us (J.P.).

Unsurprisingly, the first PCA axis (81.77% of variance) separates the avian from the non-avian Theropoda, while the second PCA axis (13.49% of variance) describes the avian variation. The interpretations that can be made from such results are however interesting.

Shape variation for PC1 (Fig. 5) ranges from a minimum very close to the non-avian Allosaurus specimen, to a maximum displaying a mix of all four birds’ attributes. This results in a relative shortening and dorsoventral compression of the preacetabular ilium, which turns from a flexed, laterally facing lateral position to a straighter, dorsally facing medial one as we move along PC1. The postacetabular ilia experience a prominent proportional lengthening, as they simultaneously widen from a dorsomedial position to a ventrolateral one. Concerning ischia, differences put forth by Procrustes superimposition reveal a relative shortening and shift from a highly medial to a fully lateral, more ventral orientation corresponding to the positive extreme of PC1. The distal pubic symphysis separates as the pubes retrovert and expand laterally alongside the ventral edge of the ilia. More generally, these results also demonstrate that non-avian theropods (i.e. Allosaurus, which is reasonably representative of most taxa) have a relatively narrow pelvis, wider at the cranial tip of the ilia, whereas birds have a relatively wider pelvis that proportionally widens toward its caudal end. Thus our methods quantitatively capture changes that have been fairly well described in qualitative anatomical descriptions (e.g. Hutchinson, 2001 and Hutchinson and Gatesy, 2000).

The hip experiences rather dramatic changes (Fig. 5 C, D, E). Acetabular cavities are found into a proportionally much more lateral position on the positive–avian–side of the axis. The acetabular opening also is proportionally much larger on the Allosaurus side of PC1 relative to the four extant birds, and this mainly involves relative displacements from the ventral halves of the iliac and ischial articular facets and from the pubes. Noticeably, along PC1 the dorsal iliac part of the acetabulum reduces mediolaterally.

The minimum of PC2 (Fig. 6) lies very close to Coturnix and the positive end of the axis lies relatively close to Troglodytes and Taeniopygia. This corresponds to a change in the shape of the preacetabular ilia from a slender, elongate outline with a broaden cranial extremity on the quail side, to a proportionally shorter, higher dorsoventrally shape with a wider base on the positive extreme of the axis. Deformation along this Principal axis shows a relative shortening of postacetabular iliac blades as they reorient in a more dorsolateral direction. Also, the ischia become slightly larger and spread relatively farther laterally on the positive extremum. Toward the positive end of the PC2 axis, the pubes relatively face more cranially, and also more laterally, though less than is observed in the ischia.

The pattern of changes occurring within the acetabulum along PC2 (Fig. 6 C, D, E) is less exaggerated than that displayed by PC1, but nonetheless may have appreciable functional importance. We observed some variation in the size of the acetabular cavity, with the PC2 minimum involving a hip with a much smaller relative diameter. There are also noteworthy changes in the orientation of the acetabular opening. On the negative end of PC2, the medial outline of the articular surface is smaller and shifted caudodorsally, resulting in a more conical opening, facing ventrolaterocranially when compared to the wider, laterally facing and more cylindrical acetabulum of the positive end. Finally, the antitrochanter also exhibits major changes. It is also less well-developed and attached only to the caudodorsal quarter of the acetabular lateral ridge on the negative end, while it is very wide and attaches over almost the entire caudodorsal half of the acetabulum on the positive side. In doing so, it ends up shifting into a relatively more dorsal orientation along PC2.

Fig. 7 represents 3D surface differences of shape along the first axis of an example of PCA realised on specimens of different geographical origin representing 38.6% of total variability. We can visualize the variable regions of the skull in 3D with a colour code. Four main examples of interpretation have to be noted.

First, the nasal part, related to olfactory sense and food foraging, is different along this axis and seems to be proportionally smaller on the positive part of the axis regarding the blue colouration of this region. This pattern may be interpreted in terms of adaptive traits of difference of olfactory capacity.

Second, a depression localised at the junction of the frontal and parietal bones distinguishes specimens along this axis.

Third, the overall form of the braincase is more or less convex and longer according to the origins of the specimens. Indeed, specimens situated on the positive part of the first axis have a proportionally longer braincase than specimens situated on the negative part which have a more globular braincase.

Fourth, the two mandibular fossa present variations that could be revealing different diets through chewing mode. In this case, these differences can be regarded in term of geometry (convexity, orientation of the articular surface) and related to the associated movement of the mandible. These last two different patterns of shape are clearly demonstrated along this axis even if further studies are needed to explain it.